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ABSTRACT 


Measurements of solar wind turbulence reveal the ubiquity of discontinuities. In 
this study, we investigate how the discontinuities, especially rotational discontinuities 
(RDs), are formed in magnetohydrodynamic (MHD) turbulence. In a simulation of 
the decaying compressive three-dimensional (3-D) MHD turbulence with an imposed 
uniform background magnetic field, we detect RDs with sharp field rotations and little 
variations of magnetic field intensity as well as mass density. At the same time, in 
the de Hoffman-Teller (HT) frame, the plasma velocity is nearly in agreement with 
the Alfvén speed, and is field-aligned on both sides of the discontinuity. We take one 
of the identified RDs to analyze in details its 3-D structure and temporal evolution. 
By checking the magnetic field and plasma parameters, we find that the identified 
RD evolves from the steepening of the Alfvén wave with moderate amplitude, and that 
steepening is caused by the nonuniformity of the Alfvén speed in the ambient turbulence. 


1. INTRODUCTION 


Since the beginning of the space age, the solar wind has been regarded as an excellent natural 
laboratory for studying the plasma turbulence. The endeavoring research over decades has revealed 
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[1994] [1996} 1996} [Horbury et al.|[2001} [Knetter et al.|/2004} [Tsurutani et al.|/2005} 
2006} [Vasquez et al.][2007} [Li][2008} [Lin et al 9000: [Sonnerup et al.|/2010} {Teh et al. 
2011; |Haaland et al.|/2012; |Malaspina & Gosling||2012; 2013} 2013). 


These discontinuities appear as large and rapid changes in properties of the plasma and magnetic 
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field, and are identified as statically advected tangential discontinuities (TDs), or propagating 
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rotational discontinuities (RDs). The TDs are characterized by small normal components of the 
magnetic field, large variations of magnetic field intensity and density jumps, while the RDs have 
large normal components of the magnetic field, but small variations of magnetic field intensity and 


density (Hudson1970). 


By measurements of the magnetic field from Pioneer 6, observed the disconti- 
nuities in the magnetic-field direction with special emphasis on their distribution in time. Based 
on interplanetary field measurements made by the vector helium magnetometers onboard Pioneer 
10 and Pioneer 11, investigated a possible dependence of the occurrence 
rate and the properties of the discontinuities on radial distance between 1 and 8.5 AU. 


(1986) surveyed the data from the Mariner 10 primary mission to study the characteris- 
tics of the discontinuities in the interplanetary magnetic field at heliographic distances of 1.0, 0.72, 


and 0.46 AU, and found an r—~!3*°-4 dependence for the daily average number of discontinuities per 
hour. With Ulysses magnetic field and plasma data obtained at radial distances ranging between 
1 and 5 AU from the Sun and at high heliographic latitudes, (Tsurutani et al. (1994) discovered 
two regions where the occurrence rate of interplanetary discontinuities is high: in stream-stream 
interaction regions and in Alfvén wave trains. To determine the normals of the discontinuities, 
explored the discontinuities measured by three spacecraft WIND, IMP 8, and 
Geotail together with the solar wind velocity measured at Geotail, and obtained quite different dis- 
tributions of the discontinuity types. With magnetic field data from the ACE spacecraft, 
extended the survey of discontinuity properties to small spread angles of the field 
vectors across the discontinuity, and found that solar wind discontinuities are far more abundant at 
small than at large spread angles. By measurements from the WIND spacecraft, 
studied the intermittent structures in solar wind turbulence, which are identified as being mostly 
rotational discontinuities (RDs) and rarely tangential discontinuities (TDs) based on the technique 
described by and (1979). carried out 
a comprehensive study of directional discontinuties and Alfvénic fluctuations in the solar wind on 
the basis of Cluster data. 


So far, there are still debates regarding the origin and nature of the discontinuities in the 
solar wind. Since RDs appear as a compressed Alfvén wave, nonlinear wave steepening has been 


suggested as the cause of its formation (Cohen & Kulsrud)1974}|/Malara & Elaoufir 1991} |Tsurutanil 
1994}|Vasquez & Hollweg}1996}}1998a)b} 2002}|Marsch & Verscharen|2011). By 


numerically calculating the evolution of an initially parallel-propagating, elliptically polarized wave 


train in a cold plasma, |Cohen & Kulsrud| (1974) were the first to investigate this possibility, and 


showed that this wave evolves to a constant-B solution with RDs that rotate the field by exactly 


180°. |Vasquez & Hollweg] (1996) continued this study, but conducted a 1.5-D hybrid numerical 


simulation study of the evolution of obliquely propagating, linearly polarized Alfvén wave trains. 
They found that large-amplitude dD / Du > 1 wave trains steepen and produce RDs which always 


rotate the field by < 180°. |Vasquez & Hollweg} (1998a)) also presented 2.5-D numerical simulations 


of a small group of nonplanar Alfvén waves to show the generation of imbedded RDs. It should 
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be noted that in their models the formation of RDs probably occurs relatively near the Sun where 
most Alfvénic fluctuations originate. 


There is another model suggesting that MHD turbulence dynamically generates these discon- 
tinuities as the solar wind flows outward. Recently, numerical simulations have been done to inves- 
tigate this assertion Ges et a. [2OT2b). 
Through analyses of MHD simulation data, examined the relationship between 


discontinuities identified by classical methods, and coherent structures identified by using inter- 
mittency statistics. They found that the two methods produce remarkably similar distributions of 
waiting times, and in fact identify many of the same events. (2009) further examined 
the link between intermittent turbulence and MHD discontinuities, directly comparing simulations 
of MHD turbulence with statistical analysis from ACE solar wind data. Their results support the 
notion that some solar wind discontinuities are consequences of intermittent turbulence. In direct 
numerical simulations of MHD turbulence with an imposed uniform magnetic field, 
investigated the statistical properties of magnetic discontinuities, and concluded that the 
discontinuities observed in the solar wind can be reproduced by MHD turbulence. However, these 
works conducted statistical studies, and did not give a clear illustration of how discontinuities, 
especially RDs, are formed in MHD turbulence. 


In the present study, we utilize a compressible 3-D MHD model to illustrate and analyze 
the forming of RDs in the turbulence. By checking the magnetic field and plasma properties, it is 
found that RD is produced by the steepening of a moderate-amplitude Alfvén wave with nonuniform 
propagating speed. The paper is organized as follows. In Section 2, a general description of the 
numerical MHD model is given. Section 3 describes the results of the numerical simulation and its 
analysis. Section 4 is reserved for the summary and discussion. 


2. NUMERICAL MHD MODEL 


The description of the plasma is given by compressible 3-D MHD, which involves a fluctu- 
ating flow velocity v(z,y,z,t), magnetic field b(z,y,z,t), density p(x,y,z,t), and temperature 
T(x,y,z,t). An uniform guide field Bo is assumed in the z—direction, so the total magnetic field 
is B= Bo +hb. 


The MHD equations are written in the following non-dimensional form: 


Op 
pe 1 
ze TN, (1) 
“ety puu + (p+ E BB| =0, (2) 


D 1 
StV: u(e +p + 5B’) - (u: B)B| = V - (B x nj), (3) 
Op 
Z +V- (uB - Bu) = nV°B , (4) 
where 1 1 
2 p 2 S 
= f { IEN = 3 
ge rh j=VxB (5) 


which corresponds to the total energy density and current density, respectively. Here, p is the mass 
density; u = (Uz, Vy, Vz) are the x, y, and z—components of velocity; p is the thermal pressure; 
B denotes the magnetic field; t is time; y = 5/3 is the adiabatic index; and 7 is the magnetic 
resistivity. 


Three independent parameters, an initial mean density pọ, a characteristic length L, and a 
characteristic plasma speed vg = 6b/.\/4mpo with 6b = (b?)!/2, are used to normalize the MHD 
equations. Other variables are normalized by their combinations. The dimensionless numbers 
appearing in the equations are the Mach number Ms = nie, where cs = \/Ypo/po is the sound 
speed, and the magnetic Reynolds number Rm = voL /n. Here, we take M, to be 0.5, consistent 
with the solar wind observations at 1 AU, and Rm to be 1000, which is limited by the available 
spatial resolution. The uniform guide field Bo is two times of the fluctuating field | db |. 


We consider periodic boundary conditions in a cube with a side length of 27 LZ and a resolution 
defined by the number of grid points which is 5123, and run a simulation from an initial state with 
kinetic and magnetic energy per unit mass (v?) = (b?) = 1. The fluctuations initially populate an 
annulus in the Fourier k—space such that 1 < k < 8, with constant amplitude and random phases 
(Dmitruk et al 20041. The initial normalized cross helicity is set to be 0.9 
such that the primordial fluctuations are highly Alfvénic. The initial density and thermal pressure 
are set to be uniform. 


To solve the equations, we employ a splitting-based finite-volume numerical scheme. The fluid 
part is solved by the Godunov-type central scheme and the magnetic part by the constrained trans- 
port approach, in conjunction with the method called second-order Monotone Upstream Schemes 
for Conservation Laws (MUSCL) for reconstruction and with the approximate Riemann solvers of 


Harten-Lax-van Leer (HLL) for calculation of the numerical fluxes (Feng et al (20111. The explicit 


second-order Runge-Kutta stepping with total variation diminishing is applied in time integration. 


3. NUMERICAL RESULTS 


Figure [I| shows the distribution of the current density in the z—direction J, in the x-z (Left) 
and z-y (Right) plane at t = 1.25. The arrows superposed on the images are the projections of 
the magnetic field vectors in the z-z (Left) and the z-y plane (Right). The black ellipses mark the 
region where the identified RD is formed. From this figure, we can see that a large-scale background 
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Fig. 1.— Distribution of the z—component of the current density J, in an z-z (Left) and z-y 


(Right) plane at t = 1.25. Superposed by arrows are the projections of the magnetic field vectors 
in the x-z (Left) and the z-y plane (Right). The black ellipses mark the region where the identified 
RD is formed. 


magnetic field is clearly present in the z—direction. As a result of the well-known anisotropic 
behavior of magnetic field fluctuations in MHD with an imposed uniform guide field 
et al.[1996), current density structures preferentially align along the guide field direction, as shown 
in the left panel, and become much more varying in the perpendicular cross section, as shown in 
the right panel. Also, J, appears to be large in magnitude at the location where the identified RD 
is formed. 


In order to detect RD, we first seek the regions with large normalized partial variance of 
increments (PVI) of the magnetic field vector. PVI in 3-D space is defined as 


| B(x + Az, y + Ay, z + Az) — B(x — Az, y — Ay, z — Az) | 
yv (| B(x + Az, y + Ay, z + Az) — Bis — Az, y — Ay, z — Az) |?) 


PVI(z,y,z) = 


where Az, Ay and Az is the grid increment in the x—, y— and z-direction, respectively. We 
first sample the magnetic field B, plasma velocity v, density p, and temperature T along a linear 
path through the region with PVI > 2, and then perform along that path the minimum variance 
analysis (MVA) using the magnetic field data to find the maximum 
variance direction (L), intermediate variance direction (M), minimum variance direction (N), and 
their corresponding eigenvalues An. A2, and A3. Finally, we check the 3-D plasma and field structure 
of the possible events. 


Figure |2| shows the magnetic field and plasma parameters along the sampling interval for an 
identified RD. In this figure, the magnetic field has been converted into the Alfvén velocity Va, 
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Fig. 2.— Magnetic field and plasma parameters along the sampling interval for an identified RD. 
Vz, Vm, and Vy denote the L, M,N components of the Alfvén velocity (black lines) and plasma 
velocity (red lines) respectively as derived from minimum variance analysis (MVA). The Horizontal 
axis displays the coordinate s along the sampling path, with s = 0 at the RD point. 


by using Va = B/,/p . We can see that the sampling series of PVI is rather close to 0, except 
near the center. The large jumps of the Alfvén velocity and plasma velocity mainly occur along 
the L direction, with Vr, jumping from -1.0 to 1.0. The N—components Vy of them are nearly 
constant, and are equal to 1.6, along with a nearly negligible jump of | Va |, the magnitude of 
Alfvén speed. Also, the Alfvén speed and plasma speed are nearly in agreement over the entire 
interval, including the jump across the RD itself. The positive correlation between them implies 
the RD propagation anti-parallel to the magnetic field. The density p, temperature T and total 
pressure P; (which consists of the summed magnetic pressure and thermal pressure) all exhibit 
relatively constant traces throughout the interval. To be noted, the jump of Vz and the large value 
of Vy, together with the relatively slight change of P, as well as p, corroborate that this event is a 
RD. 


Figure [3]exhibits the 3-D structure of the identified RD. The green lines in the left panel denote 
magnetic field lines and the yellow arrows in the right panel are plasma velocity vectors, which are 
converted into the de Hoffman-Teller (HT) frame. The light gray surface is the isosurface where 
PVI = 4, showing RD, and red, green as well as blue arrows display the x—, y—, and z—direction, 
respectively. This figure displays that in the HT frame, magnetic field and plasma velocity have 
evident normal components across the RD, and they both rotate by a certain angle. Also, the 
plasma velocity is field aligned on both sides of the discontinuity, in accordance with the Alfvénic 
nature of RD. The identified RD appears as a thin surface and its normal inclines to the z—axis. 


To see how the identified RD is formed, Figure |4| shows the isosurfaces of By (Left) and Vz 


Fig. 3.— Three-dimensional structure of the identified RD. The green lines in the left panel 
denote magnetic field lines and the yellow arrows in the right panel are plasma velocity vectors, 
which are converted into de Hoffman-Teller (HT) frame. The light gray surface is the isosurface 
where PVI = 4, showing the RD, and red, green as well as blue arrows display the r—, y—, and 
z—direction, respectively. 
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Fig. 4.— The isosurfaces of By (Left) and V; (Right) at different moments in time, with the light 
gray surface being the isosurface where PVI = 3, and exhibiting the identified RD. 
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(Right) at different moments in time. These isosurfaces are associated with the Alfvén wave as 
shown below. PVI = 3 is used here to exhibit the identified RD at these four moments, and is shown 
as the light gray surface. From this figure, it is notable that the mutual approaching of the two 
isosurfaces of By (Vz), which describes the steepening of the Alfvén wave, leads to the formation 
of a RD. At t = 0.95, the isosurface with By = —0.55 (Vy = —0.45) is relatively away from the 
isosurface with B, = 0.75 (Vs = 0.45), and the transition of By (Vz) from Bs = —0.55 (Vy = —0.45) 
to By = 0.75 (Vz = 0.45) is gentle. The isosurface where PVI = 3 is small and thin. At t = 1.10, 
the isosurfaces with D, = —0.55 (Vz = —0.45) and B, = 0.75 (Vz = 0.45) approach each other, 
and the transition between the two isosurfaces of By (Vz) becomes steep. Correspondingly, the 
isosurface where PVI = 3 grows. Afterwards, that is at t = 1.25, the two isosurfaces of By (Vz) 
are close enough to each other, and the transition between the two isosurfaces of By (Vz) becomes 
steeper. The isosurface where PVI = 3 grows larger and thicker, and the identified RD is fully 
grown. However, at t = 1.40, the isosurfaces with By = —0.55 (V = —0.45) is far away from the 
isosurface with B, = 0.75 (Vy = 0.45), and the transition between the two isosurfaces of By (Vz) 
becomes gentle again. As a result, the isosurface where PVI = 3 becomes small and thin. The 
identified RD starts to collapse. 


To understand the type of waves before the RD is formed and to investigate the process of 
the evolution of the isosurfaces mentioned above, Figure |5] exhibits the distribution of B, (the first 
row), Vz (the second row), and V4 (the third row) in the neighborhood of the identified RD in the 
plane x-z at different moments in time. Superposed by arrows are the projections of the magnetic 
field (the first row), velocity field in the HT frame (the second row), and negative Alfvén speed field 
(the third row) vectors in the z-z plane. The black ellipses mark the position where the identified 
RD is formed. Near the ellipses, the magnitudes of By and Vy are almost same (they are set to 
have the same color scales so that the same color stands for the same value), and the directions of 
the in-plane projection of the B and V vectors are almost identical. Hence in this vicinity, we have 
V = B/,/po (we recall that ou + 1), which agrees with the polarity relations of an Alfvén wave. In 
other words, the RD is detected in a neighbouring Alfvénic environment that apparently favours 
RD formation. 


Therefore, V4 can be regarded as the propagation velocity of the structures relative to the 
location where the RD forms. Hence it is significant to trace the changes of the Alfvén speed. In 
the third row in Figure |5| where the evolution of the Alfvén speed is shown, we can see that at 
t = 0.95, there is a difference of Alfvén speed across the black ellipse, which makes the layers with 
negative B, (blue in Figure, propagate faster than that with positive By (green in Figure|4). At 
t = 1.10, there is the evidence for approaching and squeezing of these two layers. The difference 
of Alfvén speed there remains. This will drive these two layers further to approach each other. 
At t = 1.25, their transition becomes sharp, and the difference of Alfvén speed nearly disappears. 
This status continues until t = 1.40, when the transition becomes gentle as a result of the faster 
propagation of the layers with positive B, than that with negative B,. It is obvious that the 
difference of Alfvén speed makes the Alfvén wave steepen, a process which forms the identified RD. 


1.5 . 1,5 
X X 


t=0.95 t=1.10 t=1.25 


Fig. 5.— Distribution of By (the first row), Vy (the second row), and V4 (the third row) in a 
subzone of the x-z plane (which passes through the identified RD) at different moments in time. 
Superposed by arrows are the projections of the magnetic field (the first row), velocity field in the 
HT frame (the second row), and negative Alfvén speed field (the third row) vectors in the x-z 
plane. The black ellipses mark the position where the identified RD is formed. 
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4. SUMMARY AND DISCUSSION 


In this study, we use a simulation of the decaying compressive 3-D MHD turbulence with an 
imposed uniform guide field as a test case to explore the formation of RDs in MHD turbulence. 
Motivated by solar wind observation at 1 AU, we consider a moderate fluctuation amplitude corre- 
sponding to 6b/Bo = 0.5 and high Alfvénic correlations with the normalized cross helicity initially 
of 0.9. A case study is thus conducted to illustrate the origin of RDs in MHD turbulence. 


The numerical simulation shows the well-known anisotropic behavior of the turbulent MHD 
field, with the current density structures preferentially aligning along the guide field direction and 
scattering in the perpendicular plane. To detect the RDs in this simulated magnetofluid, we first 
seek the regions with large PVI, then conduct a MVA by sampling the parameters along a linear 
path, and finally check the 3-D structure of the possible RD events. One clear RD is identified with 
sharp field rotations and little variations of the magnetic field intensity as well as density. At the 
same time, in the HT frame, the plasma velocity is nearly in agreement with the Alfvén speed, and 
is field aligned on both sides of the discontinuity, satisfying the Walen relation that expresses the 
Alfvénic nature of an RD. The normal direction of the identified RD inclines to the z—axis, and 
propagates anti-parallel to the guide field. 


The comprehensive information obtained by the simulation of the magnetic field and plasma 
parameters associated with the RD implies that the RD is produced by the steepening of the 
moderate-amplitude Alfvén wave with nonuniform Alfvén speed in the ambient turbulence. Before 
the RD is formed, the layers with negative By smoothly transits to the layers with positive By. 
However, there is a difference of Alfvén speed across them, which makes the layers with negative Bz 
chase after its counterpart with positive B,. As they are driven by the neigbouring turbulence to 
approach and squeeze each other, the transition between them undergoes further steepening until 
the difference of Alfvén speed nearly disappears. At the same time, the identified RD is formed. 


In this work we investigated only one RD. Certainly, there are many RDs generated in the 
simulation of decaying MHD turbulence. It will be worthwhile to conduct statistical studies to see 
whether all RDs are produced by this nonlinear steepening of an Alfvén wave. Meanwhile, TDs 
are also formed. It can be inspiring to investigate their generation mechanism in MHD turbulence. 
Furthermore, the parameters we adopted in the simulation (e.g. Mach number, cross helicity, plasma 
2) may influence the forming of RDs or TDs. In the future, we plan to conduct a parameter study 
to investigate the possible effects induced by parameter variation on the generation of discontinuties 
in MHD turbulence. 


This work is supported by NSFC under contract Nos. 41304133, 41474147, 41231069, 41222032, 
41174148, 41421003, and 41204105. The numerical calculation has been completed on computing 
system at Peking University. 
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